rm(list = ls())
#working_folder<- "/Users/bgpopescu/Library/CloudStorage/Dropbox/covid_institutions/"
working_folder<-"//Mac/Dropbox-1/covid_paper_replication/"
setwd(working_folder)

#This is to obtain:
#-table_A9
#-figure_a16a
#-figure_a16b
#-figure_a16c
#-figure_a15

library("readxl")
library("reshape")
library("stringi")
library("stringr")
library("dplyr")
library("tidyr")
library("stargazer")
library("corrplot")
library("Hmisc")
library("glue")
library("ggplot2")
library("glue")
library('gridExtra')
library('grid')
library('lfe')
library('broom')
library('gridExtra')
library("plm")
library("multiwayvcov")
library("lmtest")
library("psych")
library('lfe')
library("ggrepel")
library("lemon")
library("sandwich")
library("multiwayvcov")
library("lmtest")
library("spdep")
library("xtable")
library('gsynth')
library("pracma")
library("conleyreg")
library("vars")
library("rlist")
library("fastDummies")
library("dplyr")
library("berryFunctions")
library("mediation")
library("MatchIt")
library("tidyverse")
library("geojsonsf")


#Step1. Turning some variables to numeric
data_cov_mob = geojson_sf("./data/data_cov_mob.geojson")
data_cov_mob$pct_manufacturing<-as.numeric(data_cov_mob$pct_manufacturing)
data_cov_mob$kul<-as.numeric(data_cov_mob$kul)
data_cov_mob$nat<-as.numeric(data_cov_mob$nat)
data_cov_mob$frei<-as.numeric(data_cov_mob$frei)
data_cov_mob$soz<-as.numeric(data_cov_mob$soz)
data_cov_mob$pol<-as.numeric(data_cov_mob$pol)
data_cov_mob$inter<-as.numeric(data_cov_mob$inter)


data_cov_mob$gspo<-as.numeric(data_cov_mob$gspo)
data_cov_mob$gkul<-as.numeric(data_cov_mob$gkul)
data_cov_mob$gnat<-as.numeric(data_cov_mob$gnat)

data_cov_mob$gsoz<-as.numeric(data_cov_mob$gsoz)
data_cov_mob$gpol<-as.numeric(data_cov_mob$gpol)
data_cov_mob$ginter<-as.numeric(data_cov_mob$ginter)

#Bridging Networks
data_cov_mob$bridging_total<-data_cov_mob$gspo+ data_cov_mob$gkul+ data_cov_mob$gnat
data_cov_mob$bridging_tot_pop<-data_cov_mob$bridging_total/data_cov_mob$pop_total*1000
mean_brid<-summary(data_cov_mob$bridging_tot_pop)[4]
data_cov_mob$bridging_tot_pop_dummy<-ifelse(data_cov_mob$bridging_tot_pop>mean_brid, 1, 0)


#Bonding Networks
data_cov_mob$bonding_total<-data_cov_mob$gsoz+ data_cov_mob$gpol+ data_cov_mob$ginter
data_cov_mob$bonding_tot_pop<-data_cov_mob$bonding_total/data_cov_mob$pop_total*1000
mean_bond<-summary(data_cov_mob$bonding_tot_pop)[4]
data_cov_mob$bonding_tot_pop_dummy<-ifelse(data_cov_mob$bonding_tot_pop>mean_bond, 1, 0)

#All Networks
data_cov_mob$gvereine<-as.numeric(data_cov_mob$gvereine)
data_cov_mob$vereine_tot_pop<-data_cov_mob$gvereine/data_cov_mob$pop_total*1000


mean_ver<-summary(data_cov_mob$vereine_tot_pop)[4]
data_cov_mob$vereine_tot_pop_dummy<-ifelse(data_cov_mob$vereine_tot_pop>mean_ver, 1, 0)

data_cov_mob_week<-subset(data_cov_mob, week=="01")
data_cov_mob_week<-subset(data_cov_mob_week, select = c("GEN", 
                                                        "vereine_tot_pop",
                                                        "bridging_tot_pop", "bonding_tot_pop"))

data<-data_cov_mob
data$week_new<-as.numeric(data$week)+1
data$week_new<-paste0("0", data$week_new)  
data$week_new2<-ifelse(nchar(data$week_new)==3, str_remove(data$week_new, "^0+"), data$week_new)
data$week<-data$week_new2
data<-subset(data, select = -c(week_new, week_new2))
data$post_treat<-ifelse(as.numeric(data$week)>10, 1,0)

data$post_treat_vereine_tot_pop_dummy<-data$post_treat*data$vereine_tot_pop_dummy
data$post_treat_bridging_tot_pop_dummy<-data$post_treat*data$bridging_tot_pop_dummy
data$post_treat_bonding_tot_pop_dummy<-data$post_treat*data$bonding_tot_pop_dummy

#For Server too
# Create dummy variables:
dataf <- dummy_cols(data, select_columns = 'AGS')
dataw <- dummy_cols(dataf, select_columns = 'week')
datas <- dummy_cols(dataw, select_columns = 'SN_L')

new<-dataw %>% dplyr::select(starts_with("AGS_"))
new2<-dataw %>% dplyr::select(starts_with("week_"))
new3<-datas %>% dplyr::select(starts_with("SN_L"))


list_dummies<-names(new)
#I am removing the counties which get dropped in Stata

#Re-enable this later if necessary
list_dummies<-list_dummies[list_dummies!= "AGS_0"]
list_dummies<-list_dummies[list_dummies!= "AGS_5316"]
list_dummies<-list_dummies[list_dummies!= "AGS_8121"]
list_dummies<-list_dummies[list_dummies!= "AGS_8221"]
list_dummies<-list_dummies[list_dummies!= "AGS_8231"]
list_dummies<-list_dummies[list_dummies!= "AGS_8311"]
list_dummies<-list_dummies[list_dummies!= "AGS_8421"]
list_dummies<-list_dummies[list_dummies!= "AGS_9461"]
list_dummies<-list_dummies[list_dummies!= "AGS_13004"]
list_dummies<-list_dummies[list_dummies!= "AGS_16052"]
list_dummies<-list_dummies[list_dummies!= "AGS_16053"]
list_dummies<-list_dummies[list_dummies!= "AGS_16054"]
list_dummies<-list_dummies[list_dummies!= "AGS_16055"]
list_dummies<-list_dummies[list_dummies!= "AGS_16061"]
list_dummies<-list_dummies[list_dummies!= "AGS_16077"]

list_week_dummies<-names(new2)
list_week_dummies<-list_week_dummies[list_week_dummies!= "week_new"]
list_week_dummies<-list_week_dummies[list_week_dummies!= "week_new2"]
list_week_dummies<-list_week_dummies[list_week_dummies!= "week_numeric"]
list_week_dummies_now1<-list_week_dummies[list_week_dummies!= "week_01"]
  
#Re-enable this later if necessary
list_week_dummies2<-list_week_dummies[list_week_dummies!= "week_10"]
list_week_dummies2<-list_week_dummies2[list_week_dummies2!= "week_53"]


list_land_dummies<-names(new3)
#Re-enable this later if necessary
list_land_dummies_full<-list_land_dummies[list_land_dummies!= "SN_L"]
list_land_dummies<-list_land_dummies_full[list_land_dummies_full!= "SN_L_01"]


#Framing impacts public opinion through anxiety
#Bridging clubs impact mobility through right-wing votes.

#Step 1: We first fit the mediator model where the measure of anxiety (emo) 
#is modeled as a function of the framing treatment
#=
#Step 1: We first fit the mediator model where the measure of right-wing votes
#is modeled as a function of bridging clubs and pre-tretment covariates.

#Step 2: Next we model the outcome variable, which is a binary variable 
#indicating whether or not the participant agreed to send a letter about immigration policy 
#to his or her member of Congress.
#The explanatory variables of the outcome model include the mediator, treatment status, 
#and the same set of pre-treatment variables as those used in the mediator model.

#Step 2: Next we model the outcome variable, which is a continuous variable measuring mobility


x<-subset(datas, select = c("AGS", "bridging_tot_pop_dummy", "pct_students", "pct_pop_under_35", "uni_pop", "post_treat"))

new_data<-datas%>%
  group_by(post_treat, AGS) %>%
  dplyr::summarize(mean_mobil = mean(mean_mobil, na.rm=TRUE),
                   bridging_tot_pop_dummy = mean(bridging_tot_pop_dummy, na.rm=TRUE),
                   #pct_unempl = mean(pct_unempl, na.rm=TRUE),
                   change_afd = mean(pct_change_afd, na.rm=TRUE),
                   pct_foreign = mean(pct_foreign, na.rm=TRUE),
                   pct_turnout = mean(pct_turnout, na.rm=TRUE),
                   log_gdp_per_cap = mean(log_gdp_per_cap, na.rm=TRUE),
                   pop_den = mean(pop_den, na.rm=TRUE),
                   pct_college = mean(pct_college, na.rm=TRUE),
                   pct_pop_above_60 = mean(pct_pop_above_60, na.rm=TRUE),
                   pct_pop_under_35 = mean(pct_pop_under_35, na.rm=TRUE),
                   pct_service_sec_narrow = mean(pct_service_sec_narrow, na.rm=TRUE),
                   pct_manufacturing = mean(pct_manufacturing, na.rm=TRUE),
                   gender_ratio = mean(gender_ratio, na.rm=TRUE),
                   #pct_students = mean(pct_students, na.rm=TRUE),
                   uni_pop = mean(uni_pop, na.rm=TRUE),
                   parks_pop = mean(parks_pop, na.rm=TRUE))


#new_data$pct_pop_under_35
datax_nomis <-subset(new_data, post_treat==1)
datax_nomis <-datax_nomis[complete.cases(datax_nomis), ]



xlm1<-lm(bridging_tot_pop_dummy ~ pct_turnout, data = datax_nomis)
summary(xlm1)


xlm2<-lm(bridging_tot_pop_dummy ~ pct_turnout + log_gdp_per_cap, data = datax_nomis)
summary(xlm2)


xlm3<-lm(bridging_tot_pop_dummy ~ pct_turnout + log_gdp_per_cap + pop_den, data = datax_nomis)
summary(xlm3)


xlm4<-lm(bridging_tot_pop_dummy ~  pct_turnout + log_gdp_per_cap + pop_den + 
           pct_college, data = datax_nomis)
summary(xlm4)


xlm5<-lm(bridging_tot_pop_dummy ~ pct_turnout + log_gdp_per_cap + pop_den + 
           pct_college + pct_pop_above_60, data = datax_nomis)
summary(xlm5)


xlm6<-lm(bridging_tot_pop_dummy ~ pct_turnout + log_gdp_per_cap + pop_den + 
           pct_college + pct_pop_above_60 + pct_pop_under_35, data = datax_nomis)
summary(xlm6)


xlm7<-lm(bridging_tot_pop_dummy ~ pct_turnout + log_gdp_per_cap + pop_den + 
           pct_college + pct_pop_above_60 + pct_pop_under_35 + pct_service_sec_narrow, data = datax_nomis)
summary(xlm7)


xlm8<-lm(bridging_tot_pop_dummy ~pct_turnout + log_gdp_per_cap + pop_den + 
           pct_college + pct_pop_above_60 + pct_pop_under_35 + 
            pct_service_sec_narrow + pct_manufacturing, data = datax_nomis)
summary(xlm8)


xlm9<-lm(bridging_tot_pop_dummy ~ pct_turnout + log_gdp_per_cap + pop_den + 
            pct_college + pct_pop_above_60 + pct_pop_under_35 + 
            pct_service_sec_narrow + pct_manufacturing +
            gender_ratio, data = datax_nomis)
summary(xlm9)


xlm10<-lm(bridging_tot_pop_dummy ~ pct_turnout + log_gdp_per_cap + pop_den + 
           pct_college + 
           pct_pop_above_60 + pct_pop_under_35 +
           pct_service_sec_narrow + pct_manufacturing + 
           gender_ratio+ uni_pop, data = datax_nomis)
summary(xlm10)

xlm11<-lm(bridging_tot_pop_dummy ~ pct_turnout + log_gdp_per_cap + pop_den + 
            pct_college + 
            pct_pop_above_60 + pct_pop_under_35 +
            pct_service_sec_narrow + pct_manufacturing + 
            gender_ratio+ uni_pop + parks_pop, data = datax_nomis)
summary(xlm11)


m.out1 <- matchit(bridging_tot_pop_dummy ~ pct_turnout + log_gdp_per_cap + pop_den + 
                    pct_college + 
                    pct_pop_above_60 + pct_pop_under_35 +
                    pct_service_sec_narrow + pct_manufacturing + 
                    gender_ratio+ uni_pop + parks_pop,
                  data = datax_nomis,
                  method = "nearest", 
                  caliper = 0.9,
                  replace = F)
new<-match.data(m.out1, data = datax_nomis, distance = "prop.score")
new2<-subset(new, select = c("AGS", "prop.score"))

xlm2<-lm(bridging_tot_pop_dummy ~ pct_turnout + log_gdp_per_cap + pop_den + 
           pct_college + 
           pct_pop_above_60 + pct_pop_under_35 +
           pct_service_sec_narrow + pct_manufacturing + 
           gender_ratio+ uni_pop + parks_pop, data = new)
summary(xlm2)


multiple_res =stargazer(xlm11,
                        xlm2,
                        covariate.labels=c("Pct. Turnout",
                                           "Log GDP per Capita",
                                           "Pop. Den.",
                                           "Pct. College",
                                           "Pct. Pop. above 60",
                                           "Pct. Pop. under 35",
                                           "Pct. Empl. Hospitality and Transport",
                                           "Pct. Manufacturing",
                                           "Gender Ratio",
                                           "University-Population Ratio",
                                           "Parks-Population Ratio"),
                        keep = c("pct_turnout", 
                                 "log_gdp_per_cap", 
                                 "pop_den",
                                 "pct_college",
                                 "pct_pop_above_60",
                                 "pct_pop_under_35",
                                 "pct_service_sec_narrow",
                                 "pct_manufacturing",
                                 "gender_ratio",
                                 "uni_pop",
                                 "parks_pop"),
                        dep.var.caption = "Dependent variable: Bridging Networks",
                        dep.var.labels.include = F,
                        omit.stat=c("LL","ser","f", "rsq"), no.space=TRUE, float=FALSE)

note_latex <- "\\multicolumn{3}{l} {\\parbox[t]{8cm}{ \\footnotesize \\textit{Notes:} 
Coefficients and standard errors in parentheses from OLS regression. 
$^{*}$p$<$0.1; $^{**}$p$<$0.05; $^{***}$p$<$0.01. The unit of analysis is kreis.}}"
multiple_res [grepl("Note", multiple_res)] <- note_latex
#cat(multiple_res, file="multiple_res.tex", sep="\n")
cat(multiple_res, file="./Paper/tables/table_A9.tex", sep="\n")




#Defining empty lists
data_ame = list()
data_ade = list()
data_prop = list()
data_tot_ef = list()


for (i in unique(datas$week)){
  new_data<-subset(datas, select =c(AGS,
                                    mean_mobil,
                                    pct_change_afd,
                                    bridging_tot_pop_dummy,
                                    week,
                                    pct_turnout,
                                    log_gdp_per_cap,
                                    pop_den,
                                    pct_college,
                                    pct_pop_above_60,
                                    pct_pop_under_35,
                                    pct_service_sec_narrow,
                                    pct_manufacturing,
                                    gender_ratio,
                                    uni_pop,
                                    parks_pop))
  new_data2<-subset(new_data, week == i)
  new_data2<-left_join(new2, new_data2, by = c("AGS"="AGS"))
  new_data2 <-new_data2[complete.cases(new_data2), ]
  
  med_fit<-lm(pct_change_afd~bridging_tot_pop_dummy,
              data=new_data2)
  summary(med_fit)
  
  out_fit<-lm(mean_mobil~pct_change_afd+bridging_tot_pop_dummy,
              data=new_data2)
  summary(out_fit)

  med_out <- mediate(med_fit, out_fit, treat = "bridging_tot_pop_dummy", 
                     mediator = "pct_change_afd", boot = TRUE, sims = 1000)
  print(i)

  
  acme<-c("group", i, "acme", med_out$d0, "acme_low", med_out$d0.ci[1], "acme_high", med_out$d0.ci[2])
  ade<-c("group", i, "ade", med_out$z0, "ade_low", med_out$z0.ci[1], "ade_high", med_out$z0.ci[2])
  prop<-c("group", i, "prop", med_out$n0, "prop_low", med_out$n0.ci[1], "prop_high", med_out$n0.ci[2])
  tot_ef<-c("group", i, "tot_ef", med_out$tau.coef, "tot_low", med_out$tau.ci[1], "tot_high", med_out$tau.ci[2])
  data_ame[[i]]<-acme
  data_ade[[i]]<-ade
  data_prop[[i]]<-prop
  data_tot_ef[[i]]<-tot_ef
}

big_data1 = do.call(rbind, data_ame)
big_data2 = do.call(rbind, data_ade)
big_data3 = do.call(rbind, data_prop)
big_data4 = do.call(rbind, data_tot_ef)

big_data<-data.frame(rbind(big_data1, big_data2, big_data3, big_data4))

names(big_data)<-c("week_no_label", "week_no", "param_label", 
                   "param", "param_low_label", "param_low", "param_high_label", "param_high")
big_data$week_no<-as.numeric(big_data$week_no)                   
big_data$param<-as.numeric(big_data$param)                   
big_data$param_low<-as.numeric(big_data$param_low)                   
big_data$param_high<-as.numeric(big_data$param_high)                   
big_data<-big_data[order(big_data$param_label),]


#Step9: Adding variables for the ggplot
all_data5<-subset(big_data, param_label=="acme" | param_label == "tot_ef")

#Step9: Adding variables for the ggplot
x<-seq(-10, 42, by = 1)
x2<-c(x,x)
all_data5$observation<-NA
all_data5$observation<- x2

all_data5$observation<-as.character(all_data5$observation)
selection<-all_data5$observation[!stringr::str_detect(all_data5$observation, '-')]
all_data5$observation[!stringr::str_detect(all_data5$observation, '-')]<-paste("+", selection, sep = "")
all_data5$t<-paste("Network x \n", "(t", all_data5$observation, ")", sep = "")
all_data5$t[all_data5$t=="Network x \n(t+0)"]<-"Network x t"
all_data5$param_name<-NA
all_data5$param_name[all_data5$param_label=="acme"]<-"ACME"
all_data5$param_name[all_data5$param_label=="ade"]<-"ADE"
all_data5$param_name[all_data5$param_label=="prop"]<-"Prop. Mediated"
all_data5$param_name[all_data5$param_label=="tot_ef"]<-"Total Effect"
all_data5<-subset(all_data5, param_label=="acme" | param_label=="tot_ef")


#Step10: Making the plot
irislabs1<-seq(1,53,6)
irislabs2<-unique(all_data5[as.numeric(all_data5$week_no) %in% irislabs1,]$t)

new_mediation4<-ggplot(data = all_data5, aes(x=week_no, y=param, color = param_name)) +
  geom_point(size = 2,
             position=position_dodge(width=0.8))+
  geom_errorbar(aes(ymin = param_low, ymax = param_high),
                size = 0.3, width = 0,
                position=position_dodge(width=0.8))+
  scale_colour_manual(name = "", 
                      values=c("ACME"="black", 
                               "Total Effect"="red"))+
  scale_y_continuous(breaks = seq(-5, 12, by = 2), limits = c(-5, 12))+
  scale_x_continuous(breaks = irislabs1,
                     labels = irislabs1,
                     sec.axis = sec_axis(~.,
                                         breaks =irislabs1,
                                         labels = irislabs2))+
  geom_hline(yintercept = 0) +
  geom_vline(xintercept=10,color="red")+
  #scale_x_continuous(name = "Bandwidth (km)", breaks=seq(5,85,5)) +
  labs(x = "Week", y = "Point Estimate")+
  theme_bw() +
  theme(axis.text.x = element_text(size=12),
        axis.text.x.top  = element_text(size=10, angle = 45, margin=margin(5,5,10,5,"pt")),
        axis.text.y = element_text(size=12),
        axis.title=element_text(size=12),
        plot.title = element_text(hjust = 0.5),
        legend.position = c(1, 0),
        #Legend.position values should be between 0 and 1. c(0,0) corresponds to the "bottom left"
        #and c(1,1) corresponds to the "top right" position.
        legend.box.background = element_rect(fill='white'),
        legend.background = element_blank(),
        legend.text=element_text(size=12))
new_mediation4<-reposition_legend(new_mediation4, 'top left')
ggsave(new_mediation4, file = "paper/graphs/figure_a16a.jpg", 
       height = 15, width = 20, 
       units = "cm", dpi = 200)


#UNMATCHED

#Defining empty lists
data_ame = list()
data_ade = list()
data_prop = list()
data_tot_ef = list()


for (i in unique(datas$week)){
  new_data<-subset(datas, select =c(AGS,
                                    mean_mobil,
                                    pct_change_afd,
                                    bridging_tot_pop_dummy,
                                    week,
                                    pct_turnout,
                                    log_gdp_per_cap,
                                    pop_den,
                                    pct_college,
                                    pct_pop_above_60,
                                    pct_pop_under_35,
                                    pct_service_sec_narrow,
                                    pct_manufacturing,
                                    gender_ratio,
                                    uni_pop,
                                    parks_pop))
  new_data2<-subset(new_data, week == i)
  #new_data2<-left_join(new2, new_data2, by = c("AGS"="AGS"))
  new_data2 <-new_data2[complete.cases(new_data2), ]
  
  med_fit<-lm(pct_change_afd~bridging_tot_pop_dummy,
              data=new_data2)
  summary(med_fit)
  
  out_fit<-lm(mean_mobil~pct_change_afd+bridging_tot_pop_dummy,
              data=new_data2)
  summary(out_fit)
  
  med_out <- mediate(med_fit, out_fit, treat = "bridging_tot_pop_dummy", 
                     mediator = "pct_change_afd", boot = TRUE, sims = 1000)
  print(i)
  
  
  acme<-c("group", i, "acme", med_out$d0, "acme_low", med_out$d0.ci[1], "acme_high", med_out$d0.ci[2])
  ade<-c("group", i, "ade", med_out$z0, "ade_low", med_out$z0.ci[1], "ade_high", med_out$z0.ci[2])
  prop<-c("group", i, "prop", med_out$n0, "prop_low", med_out$n0.ci[1], "prop_high", med_out$n0.ci[2])
  tot_ef<-c("group", i, "tot_ef", med_out$tau.coef, "tot_low", med_out$tau.ci[1], "tot_high", med_out$tau.ci[2])
  data_ame[[i]]<-acme
  data_ade[[i]]<-ade
  data_prop[[i]]<-prop
  data_tot_ef[[i]]<-tot_ef
}

big_data1 = do.call(rbind, data_ame)
big_data2 = do.call(rbind, data_ade)
big_data3 = do.call(rbind, data_prop)
big_data4 = do.call(rbind, data_tot_ef)

big_data<-data.frame(rbind(big_data1, big_data2, big_data3, big_data4))

names(big_data)<-c("week_no_label", "week_no", "param_label", 
                   "param", "param_low_label", "param_low", "param_high_label", "param_high")
big_data$week_no<-as.numeric(big_data$week_no)                   
big_data$param<-as.numeric(big_data$param)                   
big_data$param_low<-as.numeric(big_data$param_low)                   
big_data$param_high<-as.numeric(big_data$param_high)                   
big_data<-big_data[order(big_data$param_label),]


#Step9: Adding variables for the ggplot

#Step9: Adding variables for the ggplot
x<-seq(-10, 42, by = 1)
x2<-c(x,x,x,x)
big_data$observation<-NA
big_data$observation<- x2

big_data$observation<-as.character(big_data$observation)
selection<-big_data$observation[!stringr::str_detect(big_data$observation, '-')]
big_data$observation[!stringr::str_detect(big_data$observation, '-')]<-paste("+", selection, sep = "")
big_data$t<-paste("Network x \n", "(t", big_data$observation, ")", sep = "")
big_data$t[big_data$t=="Network x \n(t+0)"]<-"Network x t"
big_data$param_name<-NA
big_data$param_name[big_data$param_label=="acme"]<-"ACME"
big_data$param_name[big_data$param_label=="ade"]<-"ADE"
big_data$param_name[big_data$param_label=="prop"]<-"Prop. Mediated"
big_data$param_name[big_data$param_label=="tot_ef"]<-"Total Effect"


#Step10: Making the plot
irislabs1<-seq(1,53,6)
irislabs2<-unique(big_data[as.numeric(big_data$week_no) %in% irislabs1,]$t)

all_data6<-subset(big_data, param_label=="acme" | param_label == "tot_ef")


new_mediation6<-ggplot(data = all_data6, aes(x=week_no, y=param, color = param_name)) +
  geom_point(size = 2,
             position=position_dodge(width=0.8))+
  geom_errorbar(aes(ymin = param_low, ymax = param_high),
                size = 0.3, width = 0,
                position=position_dodge(width=0.8))+
  scale_colour_manual(name = "", 
                      values=c("ACME"="black", 
                               "Total Effect"="red"))+
  scale_y_continuous(breaks = seq(-5, 12, by = 2), limits = c(-5, 12))+
  scale_x_continuous(breaks = irislabs1,
                     labels = irislabs1,
                     sec.axis = sec_axis(~.,
                                         breaks =irislabs1,
                                         labels = irislabs2))+
  geom_hline(yintercept = 0) +
  geom_vline(xintercept=10,color="red")+
  #scale_x_continuous(name = "Bandwidth (km)", breaks=seq(5,85,5)) +
  labs(x = "Week", y = "Point Estimate")+
  theme_bw() +
  theme(axis.text.x = element_text(size=12),
        axis.text.x.top  = element_text(size=10, angle = 45, margin=margin(5,5,10,5,"pt")),
        axis.text.y = element_text(size=12),
        axis.title=element_text(size=12),
        plot.title = element_text(hjust = 0.5),
        legend.position = c(1, 0),
        #Legend.position values should be between 0 and 1. c(0,0) corresponds to the "bottom left"
        #and c(1,1) corresponds to the "top right" position.
        legend.box.background = element_rect(fill='white'),
        legend.background = element_blank(),
        legend.text=element_text(size=12))
new_mediation6<-reposition_legend(new_mediation6, 'top left')
ggsave(new_mediation6, file = "paper/graphs/figure_a16b.jpg", 
       height = 15, width = 20, 
       units = "cm", dpi = 200)

ggsave(new_mediation6, file = "paper/graphs/figure_a15.jpg", 
       height = 15, width = 20, 
       units = "cm", dpi = 200)

#Making Map to show the matched counties

################################
#Making Map for turnout#
################################

Kreisgrenzen_2017 <- st_read(dsn="./data/shape_files.gdb",
                             layer="Kreisgrenzen_2017")

mapa_df<-Kreisgrenzen_2017

mapa_df$land_reg_kreis<-paste0(as.character(as.numeric(mapa_df$SN_L)), 
                               as.character(mapa_df$SN_R), 
                               as.character(as.numeric(as.character(mapa_df$SN_K))))
new2$AGS
mapa_df$AGS<-as.numeric(mapa_df$AGS)
mapa_df<-left_join(mapa_df, new2, by = c("AGS" = "AGS"))
mapa_df$sample<-NA
mapa_df$sample<-ifelse(is.na(mapa_df$prop.score), "Not Matched", "Matched")


color<-c("grey80", "grey40")

labels <-c("Matched", "Not Matched")

col = c("Matched"="grey80",
        "Not Matched"="grey40")

cols <- function(x) ggplot_build(x)$data[[1]][, "fill"]
fill_scale <- function(order) scale_fill_manual(values = col[order], 
                                                na.value = "black",
                                                name="Sample")


col[2:1]
match_map<-ggplot() + 
  #Raster works but takes some time
  #geom_raster(data=hdf,aes(X,Y,alpha=Hill)) +
  #scale_alpha(name = "Altitude in m", guide=F)  + 
  #geom_polygon(data = mapa_df, aes(x=long, y=lat), colour = "grey60", fill=NA)+
  geom_sf(data = mapa_df,
               aes(fill = sample), 
               color = "black", size = 0.4, alpha=0.8)+
  theme_bw()+
  fill_scale(2:1)+
  labs(x = "Longitude", y="Latitude")+
  #ggtitle("")+
  theme(legend.key.height=unit(0.4, "cm"),
        legend.key.width=unit(0.4, "cm"),
        legend.box.background = element_rect(fill = "white", 
                                             size=0.5, linetype="solid",
                                             colour = "black"),
        #legend.box.margin = margin(6, 6, 6, 6),
        #legend.title = element_text(colour="black", size=10, face="bold"),
        #legend.text = element_text(colour="black", size = 8),
        # Remove panel border
        panel.border = element_rect(colour = "black", fill=NA, size=0.5),  
        # Remove panel grid lines
        panel.grid.major = element_line(colour = "gray80",size=0.005),
        #panel.grid.minor = element_line(colour = "gray60",size=0.005),
        # Remove panel background
        panel.background = element_blank(),
        # Add axis line
        axis.line = element_line(colour = "grey"),
        axis.text=element_text(size=12),
        axis.title.x = element_text(size=12),
        axis.title.y = element_text(size=12),
        plot.title = element_text(hjust = 0.5, size=15),
        axis.text.x=element_text(size=12, angle=0, vjust = 0.65))+
  ggsn::scalebar(x.min = 12.5, x.max = 15,
                 #Above you mention how long the scalebar should be
                 y.min = 55, y.max = 55.5,
                 #Above you mention how thick the scalebar should be
                 dist = 100, dist_unit = "km",
                 transform = T, model = "WGS84",
                 location = "topright",
                 st.size = 4,
                 #Above you have the font size of the numbers below the scalebar
                 st.dist =0.4,
                 #Above you have the distance between the bar and the text, as a proportion of the y axis.
                 height=0.2)
match_map<-reposition_legend(match_map, 'bottom left')


ggsave(match_map, file = "paper/graphs/figure_a16c.jpg", 
       height = 20, width = 15, 
       units = "cm", dpi = 200)



